Monte Carlo Integration
   HOME

TheInfoList



OR:

In
mathematics Mathematics is an area of knowledge that includes the topics of numbers, formulas and related structures, shapes and the spaces in which they are contained, and quantities and their changes. These topics are represented in modern mathematics ...
, Monte Carlo integration is a technique for
numerical integration In analysis, numerical integration comprises a broad family of algorithms for calculating the numerical value of a definite integral, and by extension, the term is also sometimes used to describe the numerical solution of differential equations ...
using random numbers. It is a particular
Monte Carlo method Monte Carlo methods, or Monte Carlo experiments, are a broad class of computational algorithms that rely on repeated random sampling to obtain numerical results. The underlying concept is to use randomness to solve problems that might be determi ...
that numerically computes a
definite integral In mathematics, an integral assigns numbers to functions in a way that describes displacement, area, volume, and other concepts that arise by combining infinitesimal data. The process of finding integrals is called integration. Along with di ...
. While other algorithms usually evaluate the integrand at a regular grid, Monte Carlo randomly chooses points at which the integrand is evaluated. This method is particularly useful for higher-dimensional integrals. There are different methods to perform a Monte Carlo integration, such as uniform sampling,
stratified sampling In statistics, stratified sampling is a method of sampling from a population which can be partitioned into subpopulations. In statistical surveys, when subpopulations within an overall population vary, it could be advantageous to sample each s ...
,
importance sampling Importance sampling is a Monte Carlo method for evaluating properties of a particular distribution, while only having samples generated from a different distribution than the distribution of interest. Its introduction in statistics is generally att ...
, sequential Monte Carlo (also known as a particle filter), and
mean-field particle methods Mean-field particle methods are a broad class of ''interacting type'' Monte Carlo algorithms for simulating from a sequence of probability distributions satisfying a nonlinear evolution equation. These flows of probability measures can always be int ...
.


Overview

In numerical integration, methods such as the
trapezoidal rule In calculus, the trapezoidal rule (also known as the trapezoid rule or trapezium rule; see Trapezoid for more information on terminology) is a technique for approximating the definite integral. \int_a^b f(x) \, dx. The trapezoidal rule works b ...
use a deterministic approach. Monte Carlo integration, on the other hand, employs a non-deterministic approach: each realization provides a different outcome. In Monte Carlo, the final outcome is an approximation of the correct value with respective error bars, and the correct value is likely to be within those error bars. The problem Monte Carlo integration addresses is the computation of a multidimensional definite integral :I = \int_f(\overline) \, d\overline where Ω, a subset of R''m'', has volume :V = \int_d\overline The naive Monte Carlo approach is to sample points uniformly on Ω:Newman, 1999, Chap. 1. given ''N'' uniform samples, :\overline_1, \cdots, \overline_N\in \Omega, ''I'' can be approximated by : I \approx Q_N \equiv V \frac \sum_^N f(\overline_i) = V \langle f\rangle. This is because the
law of large numbers In probability theory, the law of large numbers (LLN) is a theorem that describes the result of performing the same experiment a large number of times. According to the law, the average of the results obtained from a large number of trials shou ...
ensures that : \lim_ Q_N = I. Given the estimation of ''I'' from ''QN'', the error bars of ''QN'' can be estimated by the
sample variance In probability theory and statistics, variance is the expectation of the squared deviation of a random variable from its population mean or sample mean. Variance is a measure of dispersion, meaning it is a measure of how far a set of numbers ...
using the unbiased estimate of the variance. : \mathrm(f)\equiv\sigma_N^2 = \frac \sum_^N \left (f(\overline_i) - \langle f \rangle \right )^2. which leads to : \mathrm(Q_N) = \frac \sum_^N \mathrm(f) = V^2\frac = V^2\frac. As long as the sequence : \left \ is bounded, this variance decreases asymptotically to zero as 1/''N''. The estimation of the error of ''QN'' is thus :\delta Q_N\approx\sqrt=V\frac, which decreases as \tfrac. This is
standard error of the mean The standard error (SE) of a statistic (usually an estimate of a parameter) is the standard deviation of its sampling distribution or an estimate of that standard deviation. If the statistic is the sample mean, it is called the standard error of ...
multiplied with V. This result does not depend on the number of dimensions of the integral, which is the promised advantage of Monte Carlo integration against most deterministic methods that depend exponentially on the dimension. It is important to notice that, unlike in deterministic methods, the estimate of the error is not a strict error bound; random sampling may not uncover all the important features of the integrand that can result in an underestimate of the error. While the naive Monte Carlo works for simple examples, an improvement over deterministic algorithms can only be accomplished with algorithms that use problem-specific sampling distributions. With an appropriate sample distribution it is possible to exploit the fact that almost all higher-dimensional integrands are very localized and only small subspace notably contributes to the integral. A large part of the Monte Carlo literature is dedicated in developing strategies to improve the error estimates. In particular, stratified sampling—dividing the region in sub-domains—and importance sampling—sampling from non-uniform distributions—are two examples of such techniques.


Example

A paradigmatic example of a Monte Carlo integration is the estimation of Ï€. Consider the function :H\left(x,y\right)=\begin 1 & \textx^+y^\leq1\\ 0 & \text \end and the set Ω = ˆ’1,1× ˆ’1,1with ''V'' = 4. Notice that :I_\pi = \int_\Omega H(x,y) dx dy = \pi. Thus, a crude way of calculating the value of Ï€ with Monte Carlo integration is to pick ''N'' random numbers on Ω and compute :Q_N = 4 \frac\sum_^N H(x_,y_) In the figure on the right, the relative error \tfrac is measured as a function of ''N'', confirming the \tfrac.


C example

Keep in mind that a true random number generator should be used. int i, throws = 99999, insideCircle = 0; double randX, randY, pi; srand(time(NULL)); for (i = 0; i < throws; ++i) pi = 4.0 * insideCircle / throws;


Wolfram Mathematica example

The code below describes a process of integrating the function :f(x) = \frac from 0.8 using the Monte-Carlo method in
Mathematica Wolfram Mathematica is a software system with built-in libraries for several areas of technical computing that allow machine learning, statistics, symbolic computation, data manipulation, network analysis, time series analysis, NLP, optimizat ...
: func _:= 1/(1 + Sinh *x(Log ^2); (*Sample from truncated normal distribution to speed up convergence*) Distrib _, average_, var_:= PDF ormalDistribution[average,_var_1.1*x_-_0.1.html" ;"title="verage,_var.html" ;"title="ormalDistribution[average, var">ormalDistribution[average, var 1.1*x - 0.1">verage,_var.html" ;"title="ormalDistribution[average, var">ormalDistribution[average, var 1.1*x - 0.1 n = 10; RV = RandomVariate[TruncatedDistribution[, NormalDistribution[1, 0.399, n]; Int = 1/n Total[func[RV]/Distrib[RV, 1, 0.399*Integrate[Distrib[x, 1, 0.399], ] NIntegrate[func[x], ] (*Compare with real answer*)


Recursive stratified sampling

Recursive stratified sampling is a generalization of one-dimensional
adaptive quadrature Adaptive quadrature is a numerical integration method in which the integral of a function f(x) is approximated using static quadrature rules on adaptively refined subintervals of the region of integration. Generally, adaptive algorithms are just ...
s to multi-dimensional integrals. On each recursion step the integral and the error are estimated using a plain Monte Carlo algorithm. If the error estimate is larger than the required accuracy the integration volume is divided into sub-volumes and the procedure is recursively applied to sub-volumes. The ordinary 'dividing by two' strategy does not work for multi-dimensions as the number of sub-volumes grows far too quickly to keep track. Instead one estimates along which dimension a subdivision should bring the most dividends and only subdivides the volume along this dimension. The stratified sampling algorithm concentrates the sampling points in the regions where the variance of the function is largest thus reducing the grand variance and making the sampling more effective, as shown on the illustration. The popular MISER routine implements a similar algorithm.


MISER Monte Carlo

The MISER algorithm is based on recursive
stratified sampling In statistics, stratified sampling is a method of sampling from a population which can be partitioned into subpopulations. In statistical surveys, when subpopulations within an overall population vary, it could be advantageous to sample each s ...
. This technique aims to reduce the overall integration error by concentrating integration points in the regions of highest variance.Press, 1990, pp 190-195. The idea of stratified sampling begins with the observation that for two disjoint regions ''a'' and ''b'' with Monte Carlo estimates of the integral E_a(f) and E_b(f) and variances \sigma_a^2(f) and \sigma_b^2(f), the variance Var(''f'') of the combined estimate :E(f) = \tfrac \left (E_a(f) + E_b(f) \right ) is given by, :\mathrm(f) = \frac + \frac It can be shown that this variance is minimized by distributing the points such that, :\frac = \frac Hence the smallest error estimate is obtained by allocating sample points in proportion to the standard deviation of the function in each sub-region. The MISER algorithm proceeds by bisecting the integration region along one coordinate axis to give two sub-regions at each step. The direction is chosen by examining all ''d'' possible bisections and selecting the one which will minimize the combined variance of the two sub-regions. The variance in the sub-regions is estimated by sampling with a fraction of the total number of points available to the current step. The same procedure is then repeated recursively for each of the two half-spaces from the best bisection. The remaining sample points are allocated to the sub-regions using the formula for ''Na'' and ''Nb''. This recursive allocation of integration points continues down to a user-specified depth where each sub-region is integrated using a plain Monte Carlo estimate. These individual values and their error estimates are then combined upwards to give an overall result and an estimate of its error.


Importance sampling

There are a variety of importance sampling algorithms, such as


Importance sampling algorithm

Importance sampling provides a very important tool to perform Monte-Carlo integration.Newman, 1999, Chap. 2. The main result of importance sampling to this method is that the uniform sampling of \overline is a particular case of a more generic choice, on which the samples are drawn from any distribution p(\overline). The idea is that p(\overline) can be chosen to decrease the variance of the measurement ''QN''. Consider the following example where one would like to numerically integrate a gaussian function, centered at 0, with σ = 1, from −1000 to 1000. Naturally, if the samples are drawn uniformly on the interval ˆ’1000, 1000 only a very small part of them would be significant to the integral. This can be improved by choosing a different distribution from where the samples are chosen, for instance by sampling according to a gaussian distribution centered at 0, with σ = 1. Of course the "right" choice strongly depends on the integrand. Formally, given a set of samples chosen from a distribution :p(\overline) : \qquad \overline_1, \cdots, \overline_N \in V, the estimator for ''I'' is given by : Q_N \equiv \frac \sum_^N \frac Intuitively, this says that if we pick a particular sample twice as much as other samples, we weight it half as much as the other samples. This estimator is naturally valid for uniform sampling, the case where p(\overline) is constant. The
Metropolis–Hastings algorithm In statistics and statistical physics, the Metropolis–Hastings algorithm is a Markov chain Monte Carlo (MCMC) method for obtaining a sequence of random samples from a probability distribution from which direct sampling is difficult. This seque ...
is one of the most used algorithms to generate \overline from p(\overline), thus providing an efficient way of computing integrals.


VEGAS Monte Carlo

The VEGAS algorithm approximates the exact distribution by making a number of passes over the integration region which creates the histogram of the function ''f''. Each histogram is used to define a sampling distribution for the next pass. Asymptotically this procedure converges to the desired distribution.Lepage, 1978 In order to avoid the number of histogram bins growing like ''Kd'', the probability distribution is approximated by a separable function: :g(x_1, x_2, \ldots) = g_1(x_1) g_2(x_2) \ldots so that the number of bins required is only ''Kd''. This is equivalent to locating the peaks of the function from the projections of the integrand onto the coordinate axes. The efficiency of VEGAS depends on the validity of this assumption. It is most efficient when the peaks of the integrand are well-localized. If an integrand can be rewritten in a form which is approximately separable this will increase the efficiency of integration with VEGAS. VEGAS incorporates a number of additional features, and combines both stratified sampling and importance sampling.


See also

*
Quasi-Monte Carlo method In numerical analysis, the quasi-Monte Carlo method is a method for numerical integration and solving some other problems using low-discrepancy sequences (also called quasi-random sequences or sub-random sequences). This is in contrast to the regu ...
*
Auxiliary field Monte Carlo Auxiliary-field Monte Carlo is a method that allows the calculation, by use of Monte Carlo techniques, of averages of operators in many-body quantum mechanical (Blankenbecler 1981, Ceperley 1977) or classical problems (Baeurle 2004, Baeurle 2003, ...
*
Monte Carlo method in statistical physics Monte Carlo in statistical physics refers to the application of the Monte Carlo method to problems in statistical physics, or statistical mechanics. Overview The general motivation to use the Monte Carlo method in statistical physics is to evaluat ...
*
Monte Carlo method Monte Carlo methods, or Monte Carlo experiments, are a broad class of computational algorithms that rely on repeated random sampling to obtain numerical results. The underlying concept is to use randomness to solve problems that might be determi ...
*
Variance reduction In mathematics, more specifically in the theory of Monte Carlo methods, variance reduction is a procedure used to increase the precision of the estimates obtained for a given simulation or computational effort. Every output random variable fro ...


Notes


References

* * * * * * * * * {{Cite book , last1=Robert , first1=CP , last2=Casella , first2=G , year=2004 , title=Monte Carlo Statistical Methods , publisher=Springer , edition=2nd , isbn=978-1-4419-1939-7


External links


Café math : Monte Carlo Integration
: A blog article describing Monte Carlo integration (principle, hypothesis, confidence interval)

* ttps://sites.google.com/view/chremos-group/applets/monte-carlo: Monte Carlo applet applied in statistical physics problems Monte Carlo methods Articles with example code Articles with example C code